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ABSTRACT 

Within ten nearby (d < 450 pc) Gould Belt molecular clouds we evaluate statistically the relative orientation between the magnetic field projected 
on the plane of sky, inferred from the polarized thermal emission of Galactic dust observed by Planck at 353 GHz, and the gas column density 
structures, quantified by the gradient of the column density, Vr. The selected regions, covering several degrees in size, are analysed at an effective 
angular resolution of 10' FWHM, thus sampling physical scales from 0.4 to 40 pc in the nearest cloud. The column densities in the selected regions 
range from Vr ~ 10^^ to 10^^ cm"^, and hence they correspond to the bulk of the molecular clouds. The relative orientation is evaluated pixel by 
pixel and analysed in bins of column density using the novel statistical tool called “Histogram of Relative Orientations”. Throughout this study, 
we assume that the polarized emission observed by Planck at 353 GHz is representative of the projeeted morphology of the magnetic field in 
each region, i.e., we assume a constant dust grain alignment efficiency, independent of the local environment. Within most clouds we find that the 
relative orientation changes progressively with increasing Vr, from mostly parallel or having no preferred orientation to mostly perpendicular. In 
simulations of magnetohydrodynamic turbulence in molecular clouds this trend in relative orientation is a signature of Alfvenic or sub-Alfvenic 
turbulence, implying that the magnetic field is significant for the gas dynamics at the scales probed by Planck. We compare the deduced magnetic 
field strength with estimates we obtain from other methods and discuss the implications of the Planck observations for the general picture of 
molecular cloud formation and evolution. 
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1. Introduction 

The formation and evolution of molecular clouds (MCs) and 
their substructures, from filaments to cores and eventually to 
stars, is the product of the interaction between turbulence, mag¬ 
netic fields, and gravity (Bergin & Tafalla 2007; McKee & 
Ostriker 2007). The study of the relative importance of these 
dynamical processes is limited by the observational techniques 
used to evaluate them. These limitations have been particularly 
critical when integrating magnetic fields into the general pic¬ 
ture of MC dynamics (Elmegreen & Scalo 2004; Crutcher 2012; 
Heiles & Haverkom 2012; Hennebelle & Falgarone 2012). 

There are two primary methods of measuring magnetic fields 
in the dense interstellar medium (ISM). First, observation of 
the Zeeman effect in molecular lines provides the line-of-sight 
component of the field B\\ (Crutcher 2005). Second, polariza¬ 
tion maps - in extinction from background stars and emis¬ 
sion from dust - reveal the orientation of the field averaged 
along the line of sight and projected on the plane of the sky 
(Hiltner 1949; Davis & Greenstein 1951; Hildebrand 1988; 
Planck Collaboration Int. XXI 2014). 

Analysis of the Zeeman effect observations presented 
by Crutcher et al. (2010) shows that in the diffuse ISM sam¬ 
pled by Hi lines (^h < 300 cm“^), the maximum magnetic field 
strength ^max does not scale with density. This is interpreted as 
the effect of diffuse clouds assembled by flows along magnetic 
field lines, which would increase the density but not the magnetic 
field strength. In the denser regions (^h > 300 cm“^), probed by 
OH and CN spectral lines, the same study reports a scaling of 
the maximum magnetic field strength ^ . The latter 

observation can be interpreted as the effect of isotropic contrac¬ 
tion of gas too weakly magnetized for the magnetic field to affect 
the morphology of the collapse. However, given that the obser¬ 
vations are restricted to pencil-like lines of sight and the molecu¬ 
lar tracers are not homogeneously distributed, the Zeeman effect 
measurements alone are not sufficient to determine the relative 
importance of the magnetic field at the multiple scales within 
MCs. 

The observation of starlight polarization provides an esti¬ 
mate of the projected magnetic field orientation in particular 
lines of sight. Starlight polarization observations show coher¬ 
ent magnetic fields around density structures in MCs (Pereyra 
& Magalhaes 2004; Franco et al. 2010; Sugitani et al. 2011; 
Chapman et al. 2011; Santos et al. 2014). The coherent polar¬ 
ization morphology can be interpreted as the result of dynam¬ 
ically important magnetic fields. However, these observations 
alone are not sufficient to map even the projected magnetic field 
morphology fully and in particular do not tightly constrain the 
role of magnetic fields in the formation of structure inside MCs. 

The study of magnetic field orientation within the MCs is 
possible through the observation of polarized thermal emission 
from dust. Far-infrared and submillimetre polarimetric observa¬ 
tions have been limited to small regions up to hundreds of square 
arcminutes within clouds (Fi et al. 2006; Matthews et al. 2014) 
or to large sections of the Galactic plane at a resolution of several 
degrees (Benoit et al. 2004; Bierman et al. 2011). On the scale 
of prestellar cores and cloud segments, these observations reveal 
both significant levels of polarized emission and coherent field 
morphologies (Ward-Thompson et al. 2000; Dotson et al. 2000; 
Matthews et al. 2009). 

The strength of the magnetic field projected on the plane 
of the sky can be estimated from polarization maps us¬ 
ing the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; 
Chandrasekhar & Fermi 1953). As discussed in Appendix D, it is 


assumed that the dispersion in polarization angle is entirely 
due to incompressible and isotropic turbulence. Turbulence also 
affects the motion of the gas and so broadens profiles of emission 
and absorption lines, as quantified by dispersion . In the DCF 
interpretation B^ is proportional to the ratio Application 

of the DCF method to subregions of the Taurus MC gives esti¬ 
mates of ~ 10 in low-density regions and 25 to 42 //G 
inside filamentary structures (Chapman et al. 2011). Values of 
B^ ^ 760 yuG have been found in dense parts of the Orion MC 
region (Houde et al. 2009). Because of the experimental difficul¬ 
ties involved in producing large polarization maps, a complete 
statistical study of the magnetic field variation across multiple 
scales is not yet available. 

Additional information on the effects of the magnetic field on 
the cloud structure is found by studying the magnetic field ori¬ 
entation inferred from polarization observations relative to the 
orientation of the column density structures. Patterns of rela¬ 
tive orientation have been described qualitatively in simulations 
of magnetohydrodynamic (MHD) turbulence with different de¬ 
grees of magnetization. This is quantified as half the ratio of 
the gas pressure to the mean-field magnetic pressure (Ostriker 
et al. 2001; Heitsch et al. 2001), with the resulting turbulence 
ranging from sub-Alfvenic to super-Alfvenic. Quantitative anal¬ 
ysis of simulation cubes, where the orientation of B is avail¬ 
able directly, reveals a preferred orientation relative to density 
structures that depends on the initial magnetization of the cloud 
(Hennebelle 2013; Soler et al. 2013). Using simple models of 
dust grain alignment and polarization efficiency to produce syn¬ 
thetic observations of the simulations, Soler et al. (2013) showed 
that the preferred relative orientation and its systematic depen¬ 
dence on the degree of magnetization are preserved. 

Observational studies of relative orientation have mostly re¬ 
lied on visual inspection of polarization maps ( e.g., Myers & 
Goodman 1991; Dotson 1996). This is adequate for evaluating 
general trends in the orientation of the field. However, it is lim¬ 
ited ultimately by the need to represent the field orientation with 
pseudo-vectors, because when a large polarization map is to be 
overlaid on a scalar-field map, such as intensity or column den¬ 
sity, only a selection of pseudo-vectors can be plotted. On the 
one hand, if the plotted pseudo-vectors are the result of aver¬ 
aging the Stokes parameters over a region, then the combined 
visualization illustrates different scales in the polarization and in 
the scalar field. On the other hand, if the plotted pseudo-vectors 
correspond to the polarization in a particular pixel, then the illus¬ 
trated pattern is influenced by small-scale fluctuations that might 
not be significant in evaluating any trend in relative orientation. 

Tassis et al. (2009) present a statistical study of relative ori¬ 
entation between structures in the intensity and the inferred mag¬ 
netic field from polarization measured at 350yum towards 32 
Galactic clouds in maps of a few arcminutes in size. Comparing 
the mean direction of the field to the semi-major axis of each 
cloud, they And that the field is mostly perpendicular to that 
axis. Similarly, Fi et al. (2013) compared the relative orienta¬ 
tion in 13 clouds in the Gould Belt, calculating the main cloud 
orientation from the extinction map and the mean orientation of 
the intercloud magnetic field from starlight polarization. That 
study reported a bimodal distribution of relative cloud and field 
orientations; that is, some MCs are oriented perpendicular and 
some parallel to the mean orientation of the intercloud field. In 
both studies each cloud constitutes one independent observation 


^ We use g^jJ to avoid confusion with cr^, which was introduced in 
Planck Collaboration Int. XIX (2014) as the uncertainty in the polariza¬ 
tion angle xj/. 
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Fig. 1. Magnetic field and column density measured by Planck towards the Taurus MC. The colours represent column density. The 
“drapery” pattern, produced using the line integral convolution method (LIC, Cabral & Leedom 1993), indicates the orientation of 
magnetic field lines, orthogonal to the orientation of the submillimetre polarization. 


of relative orientation, so that the statistical significance of each 
study depends on the total number of clouds observed. In a few 
regions of smaller scales, roughly a few tenths of a parsec, Koch 
et al. (2013) report a preferred orientation of the magnetic field, 
inferred from polarized dust emission, parallel to the gradient of 
the emission intensity. 

By measuring the intensity and polarization of thermal emis¬ 
sion from Galactic dust over the whole sky and down to scales 
that probe the interiors of nearby MCs, Planck^ provides an un¬ 
precedented data set from a single instrument and with a com¬ 
mon calibration scheme, for studying the morphology of the 
magnetic field in MCs and the surrounding ISM, as illustrated 
for the Taurus region in Fig. 1. We present a quantitative anal¬ 
ysis of the relative orientation in a set of nearby (J < 450 pc) 
well-known MCs to quantify the role of the magnetic field in the 
formation of density structures on physical scales ranging from 
tens of parsecs to approximately one parsec in the nearest clouds. 

The present work is an extension of previous findings, as re¬ 
ported by the Planck collaboration, on their study of the polar¬ 
ized thermal emission from Galactic dust. Previous studies in- 


^ Planck (http://www.esa.int/Planck) is a project of the 
European Space Agency (ESA) with instruments provided by two sci¬ 
entific consortia funded by ESA member states (in particular the lead 
countries Erance and Italy), with contributions from NASA (USA) and 
telescope reflectors provided by a collaboration between ESA and a sci¬ 
entific consortium led and funded by Denmark. 


elude an overview of this emission (Planck Collaboration Int. 
XIX 2014), which reported dust polarization percentages up 
to 20% at low decreasing systematically with increasing 
A^h to a low plateau for regions with A^h > lO^^cm"^. Planck 
Collaboration Int. XX (2014) presented a comparison of the po¬ 
larized thermal emission from Galactic dust with results from 
simulations of MHD turbulence, focusing on the statistics of the 
polarization fractions and angles. Synthetic observations were 
made of the simulations under the simple assumption of homo¬ 
geneous dust grain alignment efficiency. Both studies reported 
that the largest polarization fractions are reached in the most dif¬ 
fuse regions. Additionally, there is an anti-correlation between 
the polarization percentage and the dispersion of the polariza¬ 
tion angle. This anti-correlation is reproduced well by the syn¬ 
thetic observations, indicating that it is essentially caused by the 
turbulent structure of the magnetic field. 

Over most of the sky Planck Collaboration Int. XXXII 
(2014) analysed the relative orientation between density struc¬ 
tures, which is characterized by the Hessian matrix, and po¬ 
larization, revealing that most of the elongated structures (fila¬ 
ments or ridges) have counterparts in the Stokes Q and U maps. 
This implies that in these structures, the magnetic field has a 
well-defined mean direction on the scales probed by Planck. 
Furthermore, the ridges are predominantly aligned with the mag¬ 
netic field measured on the structures. This statistical trend be¬ 
comes more striking for decreasing column density and, as ex- 
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pected from the potential effects of projection, for increasing po¬ 
larization fraction. There is no alignment for the highest col¬ 
umn density ridges in the A^h ^ 10^^ cm“^ sample. Planck 
Collaboration Int. XXXIII (2014) studied the polarization prop¬ 
erties of three nearby filaments, showing by geometrical mod¬ 
elling that the magnetic field in those representative regions has 
a well-defined mean direction that is different from the field ori¬ 
entation in the surroundings. 

In the present work, we quantitatively evaluate the relative 
orientation of the magnetic field inferred from the Planck po¬ 
larization observations with respect to the gas column density 
structures, using the histogram of relative orientations (HRO, 
Soler et al. 2013). The HRO is a novel statistical tool that quan¬ 
tifies the relative orientation of each polarization measurement 
with respect to the column density gradient, making use of the 
unprecedented statistics provided by the Planck polarization ob¬ 
servations. The HRO can also be evaluated in both 3D simula¬ 
tion data cubes and synthetic observations, thereby providing a 
direct comparison between observations and the physical con¬ 
ditions included in MHD simulations. We compare the results 
of the HRO applied to the Planck observations with the results 
of the same analysis applied to synthetic observations of MHD 
simulations of super-Alfvenic, Alfvenic, and sub-Alfvenic tur¬ 
bulence. 

Thus by comparison with numerical simulations of MHD 
turbulence, the HRO provides estimates of the magnetic field 
strength without any of the assumptions involved in the DCF 
method. For comparison, we estimate using the DCF method 
and the related method described by Hildebrand et al. (2009) ( 
DCF-hSF, for DCF plus structure function) and provide a critical 
assessment of their applicability. 

This paper is organized as follows. Section 2 introduces the 
Planck 353 GHz polarization observations, the gas column den¬ 
sity maps, and the CO line observations used to derive the ve¬ 
locity information. The particular regions where we evaluate the 
relative orientation between the magnetic field and the column 
density structures are presented in Sect. 3. Section 4 describes 
the statistical tools used for the study of these relative orienta¬ 
tions. In Sect. 5 we discuss our results and their implications in 
the general picture of cloud formation. Finally, Sect. 6 summa¬ 
rizes the main results. Additional information on the selection 
of the polarization data, the estimation of uncertainties affecting 
the statistical method, and the statistical significance of the rel¬ 
ative orientation studies can be found in Appendices A, B, and 
C, respectively. Appendix D presents alternative estimates of the 
magnetic field strength in each region. 

2. Data 

2.1. Thermal dust polarization 

Over the whole sky Planck observed the linear polarization 
(Stokes Q and U) in seven frequency bands from 30 to 353 GHz 
(Planck Collaboration I 2014). In this study, we used data from 
the High Frequency Instrument (HFI, Lamarre et al. 2010) at 
353 GHz, the highest frequency band that is sensitive to polar¬ 
ization. Towards MCs the contribution of the cosmic microwave 
background (CMB) polarized emission is negligible at 353 GHz, 
making this the Planck map that is best suited to studying the 
spatial structure of the dust polarization (Planck Collaboration 
Int. XIX 2014; Planck Collaboration Int. XX 2014). 

We used the Stokes Q and U maps and the associated noise 
maps made from five independent consecutive sky surveys of 
the Planck cryogenic mission, which together correspond to the 


DR3 (delta-DXlld) internal data release. We refer to previ¬ 
ous Planck publications for the data processing, map making, 
photometric calibration, and photometric uncertainties (Planck 
Collaboration II 2014; Planck Collaboration V 2014; Planck 
Collaboration VI 2014; Planck Collaboration VIII 2014). As in 
the first Planck polarization papers, we used the International 
Astronomical Union (lAU) conventions for the polarization an¬ 
gle, measured from the local direction to the north Galactic pole 
with positive values increasing towards the east. 

The maps of Q, U, their respective variances cTq, ct^, and 
their covariance ctqu are initially at 4.'8 resolution in HEALPix 
format^ with a pixelization at Aside = 2048, which corresponds 
to an effective pixel size of 1.'7. To increase the signal-to-noise 
ratio (S/N) of extended emission, we smoothed all the maps to 
10' resolution using a Gaussian approximation to the Planck 
beam and the covariance smoothing procedures described in 
Planck Collaboration Int. XIX (2014). 

The maps of the individual regions are projected and resam¬ 
pled onto a Cartesian grid by using the gnomonic projection pro¬ 
cedure described in Paradis et al. (2012). The HRO analysis is 
performed on these projected maps. 

2.2. Column density 

We used the dust optical depth at 353 GHz (T353) as a proxy 
for the gas column density (Ah). The T353 map (Planck 
Collaboration XI 2014) was derived from the all-sky Planck in¬ 
tensity observations at 353, 545, and 857 GHz, and the IRAS ob¬ 
servations at 100/im, which were fitted using a modified black 
body spectrum. Other parameters obtained from this fit are the 
temperature and the spectral index of the dust opacity. The T353 
map, computed initially at 5' resolution, was smoothed to 10' to 
match the polarization maps. The errors resulting from smooth¬ 
ing the product T353 map, rather than the underlying data, are 
negligible compared to the uncertainties in the dust opacity and 
do not significantly affect the results of this study. 

To scale from T353 to Ah, following Planck Collaboration XI 
(2014), we adopted the dust opacity found using Galactic extinc¬ 
tion measurements of quasars, 

T353/A'H = l-2xl0-26cm^ (1) 

Variations in dust opacity are present even in the diffuse ISM 
and the opacity increases systematically by a factor of 2 from 
the diffuse to the denser ISM (Planck Collaboration XXIV 2011; 
Martin et al. 2012; Planck Collaboration XI 2014), but our re¬ 
sults do not critically depend on this calibration . 

3. Analysed regions 

The selected regions, shown in Fig. 2, correspond to nearby {d < 
450 pc) MCs, whose characteristics are well studied and can 
be used for cloud-to-cloud comparison (Poppel 1997; Reipurth 
2008). Their properties are summarized in Table 1, which in¬ 
cludes: Galactic longitude I and latitude h at the centre of the 
field; field size Al x Ab; estimate of distance; mean and maxi¬ 
mum total column densities from dust, (Ah) and max (Ah), re¬ 
spectively; and mean H 2 column density from CO. 

In the table the regions are organized from the nearest to 
the farthest in three groups: (a) regions located ai d ^ 150pc, 
namely Taurus, Ophiuchus, Lupus, Chamaeleon-Musca, and 
Corona Australis (CrA); (b) regions located at J 300 pc, 

^ Gorski et al. (2005), http: //healpix. sf.net 
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Fig. 2. Locations and sizes of the regions selected for analysis. The background map is the gas column density, derived from 
the dust optical depth at 353 GHz (Planck Collaboration XI 2014). 


Table 1. Locations and properties of the selected regions 


Region I b Al Ab Distance^ (A^h)^ Max(A^H)^ (^H 2 )^ 

n n n n [pc] [lO^icm-^J [lO^icm-^J [lO^icm-^J 


Taurus . 172.5 -14.5 15.0 15.0 140 5.4 51.9 1.6 

Ophiuchus . 354.0 17.0 13.0 13.0 140 4.4 103.3 1.1 

Lupus. 340.0 12.7 12.0 12.0 140 3.8 30.8 1.2 

Chamaeleon-Musca . 300.0 -15.0 16.0 16.0 160 2.3 29.7 1.3 

Corona Australis (CrA) 0.0 -22.0 12.0 12.0 170 1.1 40.5 1.2 


AquilaRift. 27.0 8.0 12.0 12.0 260 9.3 58.7 1.9 

Perseus. 159.0 -20.0 9.0 9.0 300 3.9 94.8 2.6 


IC5146 . 94.0 -5.5 5.0 5.0 400 3.7 22.6 1.0 

Cepheus. 110.0 15.0 16.0 16.0 440 4.2 21.3 1.2 

Orion. 212.0 -16.0 16.0 16.0 450 5.0 93.6 2.2 


^ The estimates of distances are from: Schlafiy et al. (2014) for Taurus, Ophiuchus, Perseus, IC5146, Cepheus, and Orion; Knude & Hog (1998) 
for Lupus and CrA; Whittet et al. (1997) for Chamaeleon-Musca; and Straizys et al. (2003) for Aquila Rift. 

^ Estimated from T 353 using Eq. (1) for the selected pixels defined in Appendix A. 

^ Using the line integral Wqo over -10 < V||/(km s"^) < 10 from the Dame et al. (2001) survey and Xco = L 8 x 10^^ cm“^ km"^ s. 


Aquila Rift and Perseus; and (c) regions located at J 450 pc, 

IC 5146, Cepheus, and Orion. 

Among the clouds in the first group (all shown in Fig. 3, 
left column) are Ophiuchus and Lupus, which are two regions 
with different star-forming activities but are close neighbours 
within an environment disturbed by the Sco-Cen OB associa¬ 
tion (Wilking et al. 2008; Comeron 2008). Chamaeleon-Musca 
is a region evolving in isolation, and it is relatively unperturbed 
(Luhman 2008). Taurus (see also Fig. 1), a cloud with low-mass 
star formation, appears to be formed by the material swept up by 
an ancient superbubble centred on the Cas-Tau group (Kenyon 
et al. 2008). Finally, CrA is one of the nearest regions with recent 
intermediate- and low-mass star formation, possibly formed by 


a high-velocity cloud impact on the Galactic plane (Neuhauser 
& Forbrich 2008). 

In the second group we consider Aquila Rift and Perseus, 
shown in Fig. 4 (left column). Aquila Rift is a large complex 
of dark clouds where star formation proceeds in isolated pockets 
(Eiroa et al. 2008; Prato et al. 2008). The Perseus MC is the most 
active site of on-going star formation within 300 pc of the Sun. It 
features a large velocity gradient and is located close to hot stars 
that might have impacted its structure (Bally et al. 2008). 

In the third group are IC5146, Cepheus, and Orion, shown 
in Fig. 5 (left column). IC 5146 is an MC complex in Cygnus. It 
includes an open cluster surrounded by a bright optical nebulos¬ 
ity called the Cocoon nebula, and a region of embedded lower- 
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Fig. 3. Left: Columm density map, logio(A^H/cm“^), overlaid with magnetic field pseudo-vectors whose orientations are inferred 
from the Planck 353 GHz polarization observations. The length of the pseudo-vectors is normalized so does not reflect the polariza¬ 
tion fraction. In this first group, the regions analysed are, from top to bottom, Taurus, Ophiuchus, Lupus, Chamaeleon-Musca, and 
CrA. Right: HROs for the lowest, an intermediate, and the highest Ah bin (black, blue, and red, respectively). For a given region, 
bins have equal numbers of selected pixels (see Sect. 4.1.1 and Appendix A) within the Ah ranges labelled. The intermediate bin 
corresponds to selected pixels near the blue contours in the column density images. The horizontal dashed line corresponds to the 
average per angle bin of 15 °. The widths of the shaded areas for each histogram correspond to the +1 cr uncertainties related to the 
histogram binning operation. Histograms peaking at 0° correspond to predominantly aligned with iso-An contours. Histograms 
peaking at 90° and/or -90° correspond to predominantly perpendicular to iso-A h contours. 


mass star formation known as the IC5146 Northern"^ Streamer 
(Harvey et al. 2008). The Cepheus Flare, called simply Cepheus 
in this study, is a large complex of dark clouds that seems to 
belong to an even larger expanding shell from an old supernova 
remnant (Kun et al. 2008). Orion is a dark cloud complex with 
on-going high and low mass star formation, whose structure ap¬ 
pears to be affected by multiple nearby hot stars (Bally 2008). 

When taking background/foreground emission and noise 
within these regions into account, pixels are selected for anal¬ 
ysis according to criteria for the gradient of the column density 
(Appendix A.2) and the polarization (Appendix A.2). 


^ In equatorial coordinates. 


4. Statistical study of the relative orientation of the 
magnetic field and column density structure 

4.1. Methodology 

4.1.1. Histogram of relative orientations 

We quantify the relative orientation of the magnetic field with 
respect to the column density structures using the HRO (Soler 
et al. 2013). The column density structures are characterized by 
their gradients, which are by definition perpendicular to the iso¬ 
column density curves (see calculation in Appendix B.l). The 
gradient constitutes a vector field that we compare pixel by pixel 
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Fig. 3. (continued). 
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Fig. 4. Same as Fig. 3 for the second group, Aquila Rift and Perseus. 


to the magnetic field orientation inferred from the polarization 
maps. 

In practice we use T 353 as a proxy for Ah (Sect. 2.2). The 
angle 0 between and the tangent to the T353 contours is eval¬ 
uated using^ 

0 = arctan (| V T353 x ^ |, VT353 • ^) , (2) 

where, as illustrated in Fig. 6 , V T353 is perpendicular to the tan¬ 
gent of the iso-T 353 curves, the orientation of the unit polariza¬ 
tion pseudo-vector E, perpendicular to is characterized by 
the polarization angle 

0 = ^ arctan(-f/, Q ), (3) 

and in Eq. (2), as implemented, the norm actually carries a sign 
when the range used for 0 is between -90° and 90°. 

^ In this paper we use the version of arctan with two signed arguments 
to resolve the n ambiguity in the orientation of pseudo-vectors (Planck 
Collaboration Int. XIX 2014). 


The uncertainties in 0 due to the variance of the T353 map 
and the noise properties of Stokes Q and U at each pixel are 
characterized in Appendix B. 

The gradient technique is one of multiple methods for char¬ 
acterizing the orientation of structures in a scalar field. Other 
methods, which include the Hessian matrix analysis (Molinari 
et al. 2011; Planck Collaboration Int. XXXII 2014) and the in¬ 
ertia matrix (Hennebelle 2013), are appropriate for measuring 
the orientation of ridges, i.e., the central regions of filamentary 
structures. The gradient technique is sensitive to contours and 
in that sense it is better suited to characterizing changes in the 
relative orientation in extended regions, not just on the crests of 
structures (Soler et al. 2013; Planck Collaboration Int. XXXII 
2014). Additionally, the gradient technique can sample multiple 
scales by increasing the size of the vicinity of pixels used for 
its calculation (derivative kernel; see Appendix B.l). Previous 
studies that assign an average orientation of the cloud (Tassis 
et al. 2009; Li et al. 2013) are equivalent to studying the relative 
orientation using a derivative kernel close to the cloud size. 

The selected pixels belong to the regions of each map where 
the magnitude of the gradient IVT 353 I is greater than in a dif- 






















Planck Collaboration: Probing the role of the magnetic field in the formation of structure in molecular clouds 





-90 -45 0 45 90 

<!> [deg] 


[logio(NH/cm-^)] 
00 



CD 

d 


220.6 


216.3 


212.0 
L [deg] 


207.7 


203.4 


If)" 

_ 


_ 

If) 


Orion 


z 

■D 


O) _ 



22 , 25 < log,o(NH/cm-2 
2K61< logloK/cm-i 
20,74< logi(^CNH/cm ^ 


<22.75 

<21.63 

<21.23 

^^^ 


-90 


-45 


0 

0 [deg] 


45 


90 


Fig. 5. Same as Fig. 3 for the third group, IC 5146, Cepheus, and Orion. 
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fuse reference field (Appendix A). This selection criterion aims 
at separating the structure of the cloud from the structure of 
the background using the reference field as a proxy. For each 
region the selected reference field is the region with the same 
size and Galactic latitude that has the lowest average Ar. See 
Appendix A. 1 . 

In addition to the selection on IVT 353 I, we only consider pix¬ 
els where the norms of the Stokes Q and U are larger than in 
the diffuse reference field, therefore minimizing the effect of 
background/foreground polarization external to the cloud. The 
relative orientation angle, 0 , is computed by using polarization 
measurements with a high S/N in Stokes Q and U, i.e., only con¬ 
sidering pixels with |2|/(Tg or \U\I(Tu > 3. This selection allows 
the unambiguous definition of E by constraining the uncertainty 
in the polarization angle (see Appendix A.2). 


Table 2. Fit of ^ vs. logio(AH/cm 


Region 

Chro 

Xhro 

Taurus . 

-0.53 

21.84 

Ophiuchus . 

-0.22 

22.70 

Lupus . 

-0.28 

21.72 

Chamaeleon-Musca . 

-0.51 

21.67 

Corona Australia (CrA) 

-0.11 

24.14 

Aquila Rift. 

-0.60 

22.23 

Perseus. 

-0.45 

21.76 

IC5146 . 

-0.68 

21.79 

Cepheus . 

-0.44 

21.90 

Orion. 

-0.28 

21.88 


^ See Eq. 6 and Fig. 7. 



Fig. 6. Schematic of the vectors involved in the calculation of the 
relative orientation angle 0 . 


Once we have produced a map of relative orientations for se¬ 
lected pixels following Eq. (2), we divide the map into bins of Ar 
containing an equal number of pixels and generate a histogram 
of 0 for each bin. The shape of the histogram is used to evaluate 
the preferred relative orientation in each bin directly. A concave 
histogram, peaking at 0 °, corresponds to the preferred alignment 
of with the Ar contours. A convex histogram, peaking at 
90° and/or -90°, corresponds to the preferred orientation of 
perpendicular to the Ar contours. 

The HROs in each region are computed in 25 Ar bins hav¬ 
ing equal numbers of selected pixels (10 bins in two regions with 
fewer pixels, CrA and ICS 146). The number of Ar bins is de¬ 
termined by requiring enough bins to resolve the highest Ar re¬ 
gions and at the same time maintaining enough pixels per Ar bin 
to obtain significant statistics from each histogram. The typical 
number of pixels per bin of Ar ranges from approximately 600 
in CrA to around 4 000 in Chamaeleon-Musca. We use 12 angle 
bins of width 15 °. 

The HROs of the first group of regions, the nearest dX d ^ 
150 pc, are shown in the right-hand column of Fig. 3. For the 
sake of clarity, we only present the histograms that correspond 
to three bins, namely the lowest and highest Ar and an interme¬ 
diate Ar value. The intermediate bin is the 12th (sixth in two re¬ 
gions with fewer pixels, CrA and IC 5146), and it corresponds to 
pixels near the blue contour in the image in the left-hand column 
of Fig. 3. The widths of the shaded areas for each histogram cor¬ 
respond to the 1 (T uncertainties related to the histogram binning 
operation, which are greater than the uncertainties produced by 
the variances of 2,17, and T353 (Appendix B). The sharp and nar¬ 
row features (“jitter”) in the HROs are independent of these vari¬ 
ances. They are the product of sampling the spatial correlations 
in the magnetic field over a finite region of the sky together with 


the histogram binning; these features average out when evaluat¬ 
ing the relative orientation over larger portions of the sky (Planck 
Collaboration Int. XXXII 2014). 

Although often asymmetric, most histograms reveal a 
change in the preferred relative orientation across Ar bins. The 
most significant feature in the HROs of Taurus, Ophiuchus, and 
Chamaeleon-Musca is the drastic change in relative orientation 
from parallel in the lowest Ar bin to perpendicular in the high¬ 
est Ar bin. In Lupus the behaviour at low Ar is not clear, but 
at high Ar it is clearly perpendicular. In contrast, CrA tends to 
show as parallel in the intermediate Ar bin, but no preferred 
orientation in the other Ar bins. 

The HROs of the clouds located dX d ^ 300 pc, Aquila Rift 
and Perseus, are shown in the right-hand column of Fig. 4. They 
indicate that the relative orientation is usually parallel in the low¬ 
est Ar bins and perpendicular in the highest Ar bins. 

The HROs of the third group, located 2 X d ^ 400 - 450 pc, 
IC5146, Cepheus, and Orion, are presented in the right-hand 
column of Fig. 5. In both IC 5146 and Orion the HROs for the 
highest Ar bins reveal a preferred orientation of the field perpen¬ 
dicular to the Ar contours (Orion is quite asymmetric), whereas 
the HROs corresponding to the low and intermediate Ar bins re¬ 
veal a preferred alignment of the field with Ar structures. This 
trend is also present, but less pronounced, in the Cepheus region. 


4.1.2. Histogram shape parameter^ 


The changes in the HROs are quantified using the histogram 
shape parameter defined as 


Ac Ag 
Ac + Ac 


( 4 ) 


where Ac is the area in the centre of the histogram (- 22 ?5 < 0 < 
22?5) and Ac the area in the extremes of the histogram (-90?0 < 
0 < -67? 5 and 67? 5 < 0 < 90? 0). The value of the result of 
the integration of the histogram over 45° ranges, is independent 
of the number of bins selected to represent the histogram if the 
bin widths are smaller than the integration range. 

A concave histogram corresponding to B^ mostly aligned 
with Ar contours would have ^ > 0. A convex histogram cor¬ 
responding to B^ mostly perpendicular to Ar contours would 
have ^ < 0. A flat histogram corresponding to no preferred rela¬ 
tive orientation would have ^ ~ 0 . 

The uncertainty in cr^, is obtained from 




AiAlcrl +AIctI) 

(Ae+Ae)4 


( 5 ) 
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Fig. 7. Histogram shape parameter ^ (Eqs. 4 and 5) calculated for the different A^h bins in each region. The cases ^ > 0 and ^ < 0 
correspond to the magnetic field oriented mostly parallel and perpendicular to the structure contours, respectively. For ^ ~ 0 there 
is no preferred orientation. The displayed values of Chro and Xhro were calculated from Eq. (6) and correspond to the grey dashed 
line in each plot. 


The variances of the areas, and , characterize the jitter of 

the histograms. If the jitter is large, is large compared to |^| 
and the relative orientation is indeterminate. The jitter depends 
on the number of bins in the histogram, but ^ does not. 

Figure 7 illustrates the change in ^ as a function of 
logio(A^H/cm“^) of each bin. For most of the clouds, ^ is pos¬ 
itive in the lowest and intermediate A^h bins and negative or 
close to zero in the highest bins. The most pronounced changes 


in ^ from positive to negative are seen in Taurus, Chamaeleon- 
Musca, Aquila Rift, Perseus, and IC 5146. 

The trend in ^ vs. logio(A^H/cm“^) can be fit roughly by a 
linear relation 

^ = Chro [logio(^H/cm“^) - Xhro] • (6) 

The values of Chro and Xhro in the regions analysed are sum¬ 
marized in Table 2. For the clouds with a pronounced change in 
relative orientation the slope Chro is steeper than about -0.5, 
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Fig. 8. Same as Fig. 3 for the two test regions located directly east (ChamEast) and south of Chamaeleon-Musca (ChamSouth). 


and the value Xhro for the logio(A^H/cm“^) at which the rela¬ 
tive orientation changes from parallel to perpendicular is greater 
than about 21.7. Ophiuchus, Lupus, Cepheus, and Orion are in¬ 
termediate cases, where ^ definitely does not go negative in the 
data, but seems to do so by extrapolation; these tend to have a 
shallower Chro and/or a higher Xhro- 

The least pronounced change in ^ is seen in CrA, where ^ is 
consistently positive in all bins, and the slope is very flat. We ap¬ 
plied the HRO analysis to a pair of test regions (Fig. 8) with even 
lower Ah values (logioCAn/cm"^) < 21.6; see also Fig. 11 be¬ 
low) located directly south and directly east of the Chamaeleon- 
Musca region. As in CrA, we And that is mostly parallel to 
the Ah contours, a fairly constant and no indication of pre¬ 
dominantly perpendicular relative orientation. 

4.2. Comparisons with previous studies 

The above trends in relative orientation between and the Ah 
contours in targeted MCs, where Bj^ tend to become perpendic¬ 
ular to the Ah contours at high Ah, agree with the results of the 
Hessian matrix analysis applied to Planck observations over the 


whole sky, as reported in Fig. 15 of Planck Collaboration Int. 
XXXII (2014). 

Evidence for preferential orientations in sections of some re¬ 
gions included in this study has been reported previously. The 
Taurus region has been the target of many studies. Moneti et al. 
(1984) and Chapman et al. (2011) And evidence using infrared 
polarization of background stars for a homogeneous magnetic 
fleld perpendicular to the embedded dense filamentary struc¬ 
ture. High-resolution submillimetre observations of intensity 
with Herschel And evidence of faint fllamentary structures (“stri- 
ations”), which are well correlated with the magnetic fleld orien¬ 
tation inferred from starlight polarization (Palmeirim et al. 2013) 
and perpendicular to the filament B211. Heyer et al. (2008) re¬ 
port a velocity anisotropy aligned with the magnetic fleld, which 
can be interpreted as evidence of the channeling effect of the 
magnetic fields. But the magnetic fleld in B211 and the dense 
fllamentary structures are not measured directly. As described 
above, using Planck polarization we And that B^^ is mostly 
aligned with the lowest Ah contours (20.8 < logio(AH/cm“^) < 
21.3) (see also Planck Collaboration Int. XXXIII 2014), al¬ 
though the aligned structures do not correspond to the striations. 
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which are not resolved at 10' resolution. However, we also find 
that at higher A^h the relative orientation becomes perpendicular. 

Similar studies in other regions have found evidence of stri- 
ations correlated with the starlight-inferred magnetic field ori¬ 
entation and perpendicular to the densest filamentary structures. 
The regions studied include Serpens South (Sugitani et al. 201 1), 
which is part of Aquila Rift in this study, Musca (Pereyra & 
Magalhaes 2004), and the Northern Lupus cloud (Matthews et al. 
2014). By studying Planck polarization in larger regions around 
the targets of these previous observations, we show that a sys¬ 
tematic change in relative orientation is the prevailing statistical 
trend in clouds that reach logio(Mt/cm“^) > 21.5. 


5. Discussion 

5.1. The relative orientation between and Ah structures 

5.1.1. Spatial distribution of the HRO signal 

The maps obtained using Eq. (2) characterize the relative orien¬ 
tation in each region, without assuming an organization of the 
Ah structures in ridges or filaments; HROs basically just sam¬ 
ple the orientation of Ah contours. However, the resulting maps 
of the relative orientation angle, shown for the Taurus region in 
Fig. 9, reveal that the regions that are mostly oriented parallel 
or perpendicular to the field form continuous patches, indicating 
that the HRO signal is not only coming from variations in the 
field or the Ah contours at the smallest scales in the map. 

HRO analysis on larger scales in a map can be achieved by 
considering a larger vicinity of pixels for calculating the gradi¬ 
ent VT 353 . This operation is equivalent to calculating the next- 
neighbour gradient on a map first smoothed to the scale of inter¬ 
est. The results for relative orientation after smoothing to reso¬ 
lutions of 15', 30', and 60' are illustrated for the Taurus region 
in Fig. 9. Figure 10 shows that the corresponding HROs have a 
similar behaviour for the three representative Ah values. These 
results confirm that the preferred relative orientation is not par¬ 
ticular to the smallest scales in the map, but corresponds to co¬ 
herent structures in Ah- 

A study of the preferred orientation for the whole cloud 
would be possible by smoothing the column density and polar¬ 
ization maps to a scale comparable to the cloud size. The statis¬ 
tical significance of such a study would be limited to the number 
of clouds in the sample and would not be directly comparable 
to previous studies of relative orientation of clouds, where elon¬ 
gated structures were selected to characterize the mean orienta¬ 
tion of each cloud (Li et al. 2013). 

5.1.2. Statistical significance of the HRO signal 

Our results reveal a systemic change of ^ with Ah, suggesting 
a systematic transition from magnetic field mostly parallel to 
Ah contours in the lowest Ah bins to mostly perpendicular in 
the highest Ah bins of the clouds studied. The statistical signif¬ 
icance of this change can in principle be evaluated by consid¬ 
ering the geometrical effects that influence this distribution. In 
Appendix C.l, using simulations of Q and U maps, we elimi¬ 
nate the possibility that this arises from random magnetic fields, 
random spatial correlations in the field, or the large scale struc¬ 
ture of the field. In Appendix C.2 we simply displace the Q and 
U maps spatially and repeat the analysis, showing that the sys¬ 
temic trend of ^ vs. Ah disappears for displacements greater than 

r. 


Using a set of Gaussian models, Planck Collaboration Int. 
XXXII (2014) estimated the statistical significance of this tran¬ 
sition in terms of the relative orientation between two vectors 
in 3D and their projection in 2D. As these authors emphasized, 
two vectors that are close to parallel in 3D would be projected as 
parallel in 2D for almost all viewing angles for which the pro¬ 
jections of both vectors have a non-negligible length, but on the 
other hand, the situation is more ambiguous seen in projection 
for two vectors that are perpendicular in 3D, because they can 
be projected as parallel in 2D depending on the angle of view¬ 
ing. The quantitative effects are illustrated by the simulations in 
Appendix C.3, where we consider distributions of vectors in 3D 
that are mostly parallel, mostly perpendicular, or have no pre¬ 
ferred orientation. The projection tends to make vector pairs look 
more parallel in 2D, but the distribution of relative orientations 
in 2D is quite similar though not identical to the distribution in 
3D. In particular, the signal of perpendicular orientation is not 
erased and we can conclude that two projected vectors with non- 
negligible lengths that are close to perpendicular in 2D must also 
be perpendicular in 3D. This would apply to the mostly perpen¬ 
dicular orientation for the highest bin of Ah in the Taurus region. 

5.2. Comparison with simulations of MHD turbulence 

As a complement to observations, MHD simulations can be used 
to directly probe the actual 3D orientation of the magnetic field 
B with respect to the density structures. The change in the rela¬ 
tive orientation with Ah was previously studied in MHD simula¬ 
tions using the inertia matrix and the HRO analysis (Hennebelle 
2013; Soler et al. 2013). Soler et al. (2013) showed that in 3D, 
the change in the relative orientation is related to the degree of 
magnetization. If the magnetic energy is above or comparable 
to the kinetic energy (turbulence that is sub-Alfvenic or close 
to equipartition), the less dense structures tend to be aligned 
with the magnetic field and the orientation progressively changes 
from parallel to perpendicular with increasing density. In the 
super-Alfvenic regime, where the magnetic energy is relatively 
low, there appears to be no change in relative orientation with 
increasing density, with B and density structures being mostly 
parallel. 

Soler et al. (2013) describe 2D synthetic observations of the 
MHD simulations. The synthetic observations are produced by 
integrating the simulation cubes along a direction perpendicular 
to the mean magnetic field and assuming a homogeneous dust 
grain alignment efficiency e = 1.0. The angular resolution of the 
simulation is obtained by assuming a distance d = 150 pc and 
convolving the projected map with a Gaussian beam of 10' full 
width at half maximum (FWHM). The trends in the relative ori¬ 
entation with Ah seen in 3D are also seen using the these 2D syn¬ 
thetic observations. Given that sub-Alfvenic or close to Alfvenic 
turbulence does not significantly disturb the well-ordered mean 
magnetic field, the orientation of B perpendicular to the iso¬ 
density contours is projected well for lines of sight that are not 
close to the mean magnetic field orientation. In contrast, the 
projected relative orientation produced by super-Alfvenic turbu¬ 
lence does not necessarily reflect the relative orientation in 3D 
as a result of the unorganized field structure. 

The direct comparison between the HROs of the regions 
in this study and of the synthetic observations is presented in 
Fig. 11. The trends in the relative orientation parameter, show 
that the simulation with super-Alfvenic turbulence does not un¬ 
dergo a transition in relative orientation from parallel to perpen¬ 
dicular for logio(AH/cm“^) < 23. In contrast, most of the ob¬ 
served clouds show a decrease in ^ with increasing Ah, close 
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Fig. 9. Maps of the absolute value of the relative orientation angle, |0|, in the Taurus region. These maps are produced after smoothing 
the input maps to beam FWHMs of 10', 15', 30', and 60' and then resampling the grid to sample each beam FWHM with the same 
number of pixels. The regions in red correspond to close to perpendicular to A^h structures. The regions in blue correspond to 
close to parallel to A^h structures. The black contour, corresponding to the A^h value of the intermediate contour introduced in 
Fig. 3, provides a visual reference to the cloud structure. 


to the trends seen from the simulations with Alfvenic or sub- 
Alfvenic turbulence for logio(A^H/cm“^) < 23. Furthermore, 
Xhro, the value of logio(A^H/cm“^) where ^ goes through zero, is 
near 21.7, which is consistent with the behaviour seen in the sim¬ 
ulations with super-Alfvenic or Alfvenic turbulence. Given that 
the physical conditions in the simulations (cTy = 2.0kms“^ and 
n = 500 cm“^) are typical of those in the selected regions (crvn 
is given in Table D.l), the similarities in the dependence of ^ on 
Ah suggest that the strength of the magnetic field in most of the 
regions analysed would be about the same as the mean magnetic 
fields in the Alfvenic and sub-Alfvenic turbulence simulations, 
which are 3.5 and 11/iG, respectively. However, more precise 
estimates of the magnetic field strength coming directly from 
the HROs would require further sampling of the magnetization 
in the MHD simulations and detailed modelling of the effects of 
the line-of-sight integration. 


Indirectly, the presence of a Ah threshold in the switch in 
preferred relative orientation between and the Ah structures 
hints that gravity plays a significant role. By contrast, for the 
regions with low average Ah (i.e., CrA and the two test regions; 
Fig. 11), there is little change in ^ and certainly no switch in the 
preferred relative orientation to perpendicular. 

5.3. Physics of the relative orientation 

The finding of dense structures mostly perpendicular to the 
magnetic field (and the small mass-to-fiux ratios discussed in 
Appendix D.4) suggests that the magnetic field in most of the 
observed regions is significant for the structure and dynamics. 
However, discerning the underlying geometry is not obvious. 
As one guide, a frozen-in and strong interstellar magnetic field 
would naturally cause a self-gravitating, static cloud to become 
oblate, with its major axis perpendicular to the field lines, be- 
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Fig. 10. HROs of the Taurus region after smoothing the input maps to beam FWHMs of 15', 30', and 60', shown from left to right, 
respectively. 


cause gravitational collapse would be restricted to occurring 
along field lines (Mouschovias 1976a,b). In the case of less 
dense structures that are not self-gravitating, the velocity shear 
can stretch matter and field lines in the same direction, thereby 
producing aligned structures, as discussed in Hennebelle (2013) 
and Planck Collaboration Int. XXXII (2014). 

If the MCs are isolated entities and the magnetic field is 
strong enough to set a preferred direction for the gravitational 
collapse, the condensations embedded in the cloud are not very 
likely to have higher column densities than their surroundings 
(Nakano 1998). This means that the formation of dense sub¬ 
structures, such as prestellar cores and stars, by gravitational 
collapse would be possible only if the matter decouples from the 
magnetic field. This is possible through the decoupling between 
neutral and ionized species (ambipolar diffusion, Mouschovias 
1991; Li & Houde 2008) or through removal of magnetic flux 
from clouds via turbulent reconnection (Lazarian & Vishniac 
1999; Santos-Lima et al. 2012). 

Alternatively, if we regard MCs not as isolated entities but 
as the result of an accumulation of gas by large-scale flows 
(Ballesteros-Paredes et al. 1999; Hartmann et al. 2001; Koyama 
& Inutsuka 2002; Audit & Hennebelle 2005; Heitsch et al. 
2006), the material swept up by colliding flows may eventually 
form a self-gravitating cloud. If the magnetic field is strong the 
accumulation of material is favoured along the magnetic field 
lines, thus producing dense structures that are mostly perpendic¬ 
ular to the magnetic field. The inflow of material might even¬ 
tually increase the gravitational energy in parts of the cloud, 
thereby producing supercritical structures such as prestellar 
cores. 

For supersonic turbulence in the ISM and MCs, density 
structures can be formed by gas compression in shocks. If 
the turbulence is strong with respect to the magnetic field 
(super-Alfvenic), gas compression by shocks is approximately 
isotropic; because magnetic flux is frozen into matter, field lines 
are dragged along with the gas, forming structures that tend 
to be aligned with the field. If the turbulence is weak with re¬ 
spect to the magnetic field (sub-Alfvenic), the fields produce a 
clear anisotropy in MHD turbulence (Sridhar & Goldreich 1994; 
Goldreich & Sridhar 1995; Matthaeus et al. 2008; Banerjee et al. 
2009) and compression by shocks that is favored to occur along 
the magnetic field lines, creating structures perpendicular to the 
field. The cold phase gas that constitutes the cloud receives no 
information about the original flow direction because the mag¬ 


netic field redistributes the kinetic energy of the inflows (Heitsch 
et al. 2009; Inoue & Inutsuka 2009; Burkhart et al. 2014). This 
seems to be the case in most of the observed regions, where the 
mostly perpendicular relative orientation between the magnetic 
field and the high column density structures is an indication of 
the anisotropy produced by the field. 

The threshold of logio(Mi/cm“^) 21.7 above which the 
preferential orientation of switches to being perpendicular 
to the Ah contours is intriguing. Is there a universal threshold 
column density that is independent of the particular MC envi¬ 
ronment and relevant in the context of star formation? In princi¬ 
ple, this threshold might be related to the column density of fila¬ 
ments at which substructure forms, as reported in an analysis of 
Herschel observations (Arzoumanian et al. 2013), but the Planck 
polarization observations leading to do not fully resolve such 
filamentary structures. In principle, this threshold might also be 
related to the column density at which the magnetic field starts 
scaling with density, according to the Zeeman effect observa¬ 
tions of B\\ (Figure 7 in Crutcher 2012). However, establishing 
such relationships requires further studies with MHD simula¬ 
tions to identify what densities and scales influence the change in 
relative orientation between B^ and Ah structures and to model 
the potential imprint in B\\ observations and in observations 
to be carried out at higher resolution. 

5.4. Effect of dust grain alignment 

Throughout this study we assume that the polarized emission ob¬ 
served by Planck at 353 GHz is representative of the projected 
morphology of the magnetic field in each region; i.e., we as¬ 
sume a constant dust grain alignment efficiency (e) that is in¬ 
dependent of the local environment. Indeed, observations and 
MHD simulations under this assumption (Planck Collaboration 
Int. XIX 2014; Planck Collaboration Int. XX 2014) indicate that 
depolarization effects at large and intermediate scales in MCs 
might arise from the random component of the magnetic field 
along the line of sight. On the other hand, the sharp drop in 
the polarization fraction at Ah > 10^^ cm“^ (reported in Planck 
Collaboration Int. XIX 2014), when seen at small scales, might 
be interpreted in terms of a decrease of e with increasing column 
density (Matthews et al. 2001; Whittet et al. 2008). 

A leading theory for the process of dust grain alignment in¬ 
volves radiative torques by the incident radiation (Lazarian & 
Hoang 2007; Hoang & Lazarian 2009; Andersson 2015). A crit- 
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Fig. 11. Histogram shape parameter ^ (Eqs. 4 and 5) calculated 
for the different A^h bins in each region. Top: relative orienta¬ 
tion in synthetic observations of simulations with super-Alfvenic 
(blue), Alfvenic (green), and sub-Alfvenic (red) turbulence, as 
detailed in Soler et al. (2013). Middle: relative orientation in 
the regions selected from the Planck all-sky observations, from 
Fig. 7. The blue data points correspond to the lowest Ah regions 
(CrA and the test regions in Fig. 8, ChamSouth and ChamFast) 
and the orange correspond to the rest of the clouds. Bottom: 
comparison between the trends in the synthetic observations (in 
colours) and the regions studied (grey). The observed smooth 
transition from mostly parallel (^ > 0) to perpendicular (^ < 0) 
is similar to the transition in the simulations for which the turbu¬ 
lence is Alfvenic or sub-Alfveic. 


ical parameter for this mechanism is the ratio between the dust 
grain size and the radiation wavelength. As the dust column den¬ 
sity increases, only the longer wavelength radiation penetrates 
the cloud and the alignment decreases. Grains within a cloud 
(without embedded sources) should have lower e than those at 
the periphery of the same cloud. There is evidence for this from 
near-infrared interstellar polarization and submillimetre polar¬ 
ization along lines of sight through starless cores (Jones et al. 
2015), albeit on smaller scales and higher column densities than 
considered here. If e inside the cloud is very low, the observed 


polarized intensity would arise from the dust in the outer layers, 
tracing the magnetic field in the “skin” of the cloud. Then the 
observed orientation of is not necessarily correlated with the 
column density structure, which is seen in total intensity, or with 
the magnetic field deep in the cloud. 

Soler et al. (2013) presented the results of HRO analysis on 
a series of synthetic observations produced using models of how 
6 might decrease with increasing density. They showed that with 
a steep decrease there is no visible correlation between the in¬ 
ferred magnetic field orientation and the high-An structure, cor¬ 
responding to nearly flat HROs. 

The HRO analysis of MCs carried out here reveals a correla¬ 
tion between the polarization orientation and the column density 
structure. This suggests that the dust polarized emission samples 
the magnetic field structure homogeneously on the scales being 
probed at the resolution of the Planck observations or, alterna¬ 
tively, that the field deep within high-A h structures has the same 
orientation of the field in the skin. 

6. Conclusions 

We have presented a study of the relative orientation of the mag¬ 
netic field projected on the plane of the sky as inferred 

from the Planck dust polarized thermal emission, with respect to 
structures detected in gas column density (Ah). The relative ori¬ 
entation study was performed by using the histogram of relative 
orientations (HRO), a novel statistical tool for characterizing ex¬ 
tended polarization maps. With the unprecedented statistics of 
polarization observations in extended maps obtained by Planck, 
we analyze the HRO in regions with different column densities 
within ten nearby molecular clouds (MCs) and two test fields. 

In most of the regions analysed we find that the relative ori¬ 
entation between B^ and Ah structures changes systematically 
with Ah from being parallel in the lowest column density ar¬ 
eas to perpendicular in the highest column density areas. The 
switch occurs at logio(AH/cm“^) 21.7. This change in relative 
orientation is particularly significant given that projection tends 
to produce more parallel pseudo-vectors in 2D (the domain of 
observations) than exist in 3D. 

The HROs in these MCs reveal that most of the high Ah 
structures in each cloud are mostly oriented perpendicular to the 
magnetic field, suggesting that they may have formed by mate¬ 
rial accumulation and gravitational collapse along the magnetic 
field lines. According to a similar study where the same method 
was applied to MHD simulations, this trend is only possible if 
the turbulence is Alfvenic or sub-Alfvenic. This implies that the 
magnetic field is significant for the gas dynamics on the scales 
sampled by Planck. The estimated mean magnetic field strength 
is about 4 and 12yuG for the case of Alfvenic and sub-Alfvenic 
turbulence, respectively. 

We also estimate the magnetic field strength in the MCs stud¬ 
ied using the DCF and DCF-rSF methods. The estimates found 
seem consistent with the above values from the HRO analysis, 
but given the assumptions and systematic effects involved, we 
recommend that these rough estimates be treated with caution. 
According to these estimates the analysed regions appear to be 
magnetically sub-critical. This result is also consistent with the 
conclusions of the HRO analysis. Specific tools, such as the DCF 
and DCF-hSF methods, are best suited to the scales and physical 
conditions in which their underlying assumptions are valid. The 
study of large polarization maps covering multiple scales calls 
for generic statistical tools, such as the HRO, for characterizing 
their properties and establishing a direct relation to the physical 
conditions included in MHD simulations. 
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The study of the structure on smaller scales is beyond the 
scope of this work, however, the presence of gravitationally 
bound structures within the MCs, such as prestellar cores and 
stars, suggests that the role of magnetic fields is changing on 
different scales. Even if the magnetic field is important in the 
accumulation of matter that leads to the formation of the cloud, 
effects such as matter decoupling from the magnetic field and 
the inflow of matter from the cloud environment lead to the for¬ 
mation of magnetically supercritical structures on smaller scales. 
Further studies will help to identify the dynamical processes that 
connect the MC structure with the process of star formation. 
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Appendix A: Selection of data 

The HRO analysis is applied to each MC using common criteria 
for selecting the areas in which the relative orientation is to be 
assessed. 

A.1. Gradient mask 

The dust optical depth, T353, observed in each region, can be in¬ 
terpreted as 

^353 “ ^53 ^353 ’ (A.l) 

where is the optical depth of the MC, T333 is the optical 
depth of the diffuse regions behind and/or in front of the cloud 
(background/foreground), and ^^-353 the noise in the optical depth 

map with variance . 

The gradient of the optical depth can be then written as 

Vr?r3 = VT^53 + V(r^«3+,5,3,3). (A.2) 

We quantify the contribution of the background/foreground and 
the noise, V(t 333 + ^t 353 )5 by evaluating VT 353 in a reference 
field with lower submillimetre emission. Given that the domi¬ 
nant contribution to the background/foreground gradient would 
come from the gradient in emission from the Galactic plane, for 
each of the regions analysed we chose a reference field of the 
same size at the same Galactic latitude and with the lowest aver¬ 
age A^h in the corresponding latitude band. We compute the av¬ 
erage of the gradient norm in the reference field, and 

use this value as a threshold for selecting the regions of the map 
where VT353 carries significant information about the structure 
of the cloud. We note that this threshold includes a contribution 
from the noise . The HROs presented in this study corre¬ 
spond to regions in each field where IVT 353 I > (iVr^^l^. 

A.2. Polarization mask 

The total Stokes parameters Q and U measured in each region 
can be interpreted as 

gOBS ^ gMC qBG ^ ^OBS ^ jjMC jjBG ^ ^ 


where and correspond to the polarized emission from 
the MC, and correspond to the polarized emission from 
the diffuse background/foreground, and 6q and 6u are the noise 
contributions to the observations, such that the variances = 

S^Q and (T^ = 6^jj. 

As in the treatment of the gradient, we estimate the contri¬ 
butions of the background/foreground polarized emission and 
the noise using the rms of the Stokes parameters in the same 
reference field, and "^be HROs presented in this 

study correspond to pixels in each region where |g| > 2 g^^g 
or \U\ > 2U^^. The “or” conditional avoids biasing the se¬ 
lected values of polarization. This first selection criterion pro¬ 
vides a similar sample to the alternative coordinate-independent 
criterion > 2 V(2ms)^ + This first criterion 

aims to distinguish between the polarized emission coming from 
the cloud and the polarized emission coming from the back¬ 
ground/foreground estimated in the reference regions. 

Additionally, as a second criterion, our sample is restricted 
to polarization measurements where |g| > Sctq or \U\ > 3cru. 
This aims to select pixels where the uncertainty in the polar¬ 
ization angle is smaller than the size of the angle bins used for 
the constructions of the HRO (Serkowski 1958; Montier et al. 
2015). In terms of the total polarized intensity, P = 
and following equations. B.4 and B.5 in Planck Collaboration 
Int. XIX (2014), the second criterion corresponds to P/crp > 3 
and uncertainties in the polarized orientation angle cr^ < 10 °. 

The fractions of pixels considered in each region, after apply¬ 
ing the selection criteria described above, are summarized in 
Table A.l. The largest masked portions of the regions corre¬ 
spond to the gradient mask, which selects mostly those areas 
of column density above the mean column density of the back¬ 
ground/foreground The polarization mask provides an in¬ 

dependent criterion that is less restrictive. The intersection of 
these two masks selects the fraction of pixels considered for the 
HRO analysis. 

Appendix B: Construction of the histogram of 
relative orientations and related uncertainties 

The HROs were calculated for 25 column density bins having 
equal numbers of selected pixels (10 bins in two regions with 
fewer pixels, CrA and IC 5146). For each of these HROs we use 
12 angle bins of width 15 ° (see Sect. 4.1.1). 

8.1. Calculation of VT 353 and the uncertainty of its orientation 

The optical depth gradient (VT 353 ) is calculated by convolving 
the T 353 map with a Gaussian derivative kernel (Soler et al. 
2013), such that VT 353 corresponds to 

VT 353 = (Gjc 0 T 353 )/ -r (Gy 0 T 353)7 = gxi-^ gyj, (B.l) 

where Gx and Gy are the kernels calculated using the x- and y- 
derivatives of a symmetric two-dimensional Gaussian function. 
The orientation of the iso-T 353 contour is calculated from the 
components of the gradient vector, 

0 = arctan , gyj . (B.2) 

Because the calculation of the gradient through convolution 
is a linear operation, the associated uncertainties can be calcu¬ 
lated using the same operation,so that 

^gx = ^T353 5 ^^3, = Gy <S) ^r353 , (B.3) 
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Table A.l. Selection of data^. 


Region 

(iv-) 

/v 

/pol 

/tot 


[lO^Ocm-^] 

[%] 

[%] 

[%] 

Taurus . 

6.2 

66 

78 

28 

Ophiuchus . 

5.6 

65 

82 

31 

Lupus . 

9.6 

65 

67 

24 

Chamaeleon-Musca . 

6.3 

49 

83 

31 

Corona Australia (CrA) 

5.4 

34 

40 

4 

AquilaRift. 

18.4 

48 

96 

32 

Perseus. 

5.3 

60 

59 

16 

IC5146 . 

26.4 

38 

93 

29 

Cepheus . 

7.5 

66 

80 

36 

Orion. 

6.1 

63 

67 

24 

ChamEast. 

6.3 

33 

38 

10 

ChamSouth. 

4.9 

36 

42 

13 


^ Mean column density of the background/foreground for each region 
estimated from a reference field at the same Galactic latitude; 

percentage fy of all pixels where IVT 353 I > percentage /poi 

of all pixels where \Q\ > or \U\ > (first polarization 

criterion) and where \Q\ > Sctq or \U\ > 3cru (second polarization 
criterion); and percentage /tot of all pixels used for the HRO analysis. 


from which we obtain the crl = 61 and cri = 6l . 

6x 6x 6y sy 

The standard deviation of the angle 6 can be written as 


ere = 



which corresponds to 


0-0 




4 




+ gla-^ 


8y 


(B.4) 


(B.5) 


In the application discussed here, the standard deviations in the 
T 353 map within the selected areas are much less than a few per¬ 
centage points, so their effect on the estimate of the orientation 
of the gradient is negligible. 


B.2. Uncertainties affecting the characterization of relative 
orientations within MCs 

B.2.1. Uncertainties in the construction of the histogram 

To estimate the uncertainty associated with the noise in Stokes 
Q and U, we produce 1000 noise realizations, Qr and Ur, using 
Monte Carlo sampling. We assume that the errors are normally 
distributed and are centred on the measured values Q and U with 
dispersions ctq and (Tu (Planck Collaboration II 2014; Planck 
Collaboration VI 2014; Planck Collaboration V 2014; Planck 
Collaboration VIII 2014). Given that ctqu is smaller than cTq 

and (T^, it is justified to generate Qr and Ur independently of 
each other. We then introduce Qr and Ur in the analysis pipeline 
and compute the HRO using the corresponding T353 map in each 
region. The results, presented in Fig. B.l for the Taurus region, 
show that the noise in Q and U does not critically affect the 
shape of the HROs or the trend in Together, the low noise in 
the maps of T353 and the selection criteria for the polarization 
measurements ensure that the HRO is well determined. 



0 [deg] 


Fig. B.l. HROs in the Taurus region that correspond to the indi¬ 
cated Nu bins. The plotted values are obtained using the original 
T353 map at 10 ' resolution and maps of the Stokes parameters 
Qr and Ur, which correspond to 1000 random noise realizations. 
Each realization is generated using a Gaussian probability den¬ 
sity function centred on the measured values Q and U with vari¬ 
ances cTq and cr^. 


Another source of uncertainty in the HRO resides in the his¬ 
togram binning process. The variance in the ^th histogram bin is 
given by 


where hk is the number of samples in the ^th bin, and /ttot is the 
total number of samples. 

Of the two independent sources of uncertainty above, we find 
that the largest contribution comes from the binning process, so 
that these are the ones shown as the shaded uncertainty ranges in 
all figures of HROs, for example for Taurus in Fig. 3. Because 
of the large number of samples in each histogram bin, the uncer¬ 
tainties in the HRO do not significantly affect the results of this 
study. 


B.2.2. Uncertainties in the histogram shape parameter ^ 

As in the case of the HRO, the uncertainty in as defined in 
Eq. (4), can be quantified using the random realizations intro¬ 
duced in the previous section. Figure B.2 shows the dependence 
of ^ on log 10 Ah obtained using Qr and Ur for the Taurus re¬ 
gion. The small variations around the trend line indicate that the 
uncertainties in Q and U do not significantly affect the trends 
discussed in this study, as expected from the behaviour of the 
histograms presented in Fig. B.l. The main source of uncertainty 
in the estimation of ^ is related to the histogram binning, char¬ 
acterized by the error bars calculated using Eq. (5). As seen in 
Fig. 7 and reproduced in Fig. C.l, these are much larger than the 
dispersion of the values of ^ in Fig. B.2. 
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Fig.B.2. Histogram shape parameter, as a function of 
logio(A^H/cm“^) in the Taurus region. The values are obtained 
using the T 353 map at 10' resolution and maps of the Stokes pa¬ 
rameters Qr and Ur that correspond to 1000 random-noise re¬ 
alizations. Each realization is generated using a Gaussian prob¬ 
ability density function centred on the measured values Q and 
f/with variances and cr^. By joining the values at each A^h 
bin, we find a trend very close to the black line in Fig. 7, with 
little dispersion from the noise; much larger are the uncertainties 
in evaluating ^ at each A^h bin, as given in Fig. 7 but not shown 
here. 


a Gaussian distribution and unit-length polarization pseudo¬ 
vectors. The polarization angle distribution is centred at if/o = 0° 
with a standard deviation = 45°. The results of this numerical 
experiment are shown in Fig. C.2 for the Taurus region. 



Fig. C.2. Like Fig. C.l but for values of ^ (red) obtained using 
1000 random realizations of Q and U maps corresponding to a 
Gaussian distribution of 0 ^ centred on 0o = 0° and with standard 
deviation = 45°. 


Appendix C: Statistical significance of the HRO 
signal 

C.1. A product of chance? 

To investigate various potential sources of the signal found in 
the HROs we use the T353 map at 10 ' resolution in each region in 
combination with Q and U maps that are produced with different 
recipes, each with 1000 realizations. 

To eliminate random fields, we use Q and U maps produced 
with a random realization of 0 with a uniform distribution and 
with unit-length polarization pseudo-vectors. The results of this 
numerical experiment are shown in Fig. C.l for the Taurus re¬ 
gion. For each of the 1000 realizations, we join the values at 
each A^h bin to show the trend lines in order to compare with the 
lines in Fig. 7. 


To eliminate random spatial correlations, we use Q and U 
maps produced from random realizations with a power spectrum 
P(k) oc For the spectral indices we adopted, = -1.5, 
-2.5, and -3.5, correlations are introduced in the orientation of 
the polarization pseudo-vectors that are independent of the struc¬ 
ture of matter. This test evaluates the statistical significance of 
the random sampling of spatial correlations in the magnetic field 
implied when we calculate the HROs in a finite region of the sky. 
The results of this numerical experiment are shown in Fig. C.3 
for the Taurus region. 

The results of these numerical experiments show that the 
trends with A^h found in the HROs and ^ do not arise by chance 
from these potential sources. 

C.2. A product of random correlations in the polarization 
maps? 


q 
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Fig. C.l. Histogram shape parameter, as a function of 
logio(A^H/cm“^) in the Taurus region. Blue lines join the ^ val¬ 
ues obtained using the T353 map at 10 ' resolution and each of 
1000 random realizations of Q and U maps corresponding to a 
uniform distribution of 0 . Results in black and grey are from the 
analysis of the Planck data, as reported in Fig. 7. 



21.0 21.5 22.0 

login (NH/cm"*) 


22.5 


To eliminate the large-scale magnetic field as the source we 
use Q and U maps produced with random realizations of ij/ with 


The random realizations of the magnetic field presented above 
are useful for characterizing the behaviour of However, they 
are bound to produce very different statistics compared to those 
in the real observations. In reality the orientations of have 
some non-trivial correlation structure in the map. In principle, 
the difference in the HROs might be due to the different correla¬ 
tions of the observed and of the random realizations of B^. 

To evaluate whether the signal present in the HROs arises 
from an actual physical relationship between B^ and T353, we 
introduce randomness by shifting the Stokes Q and U maps with 
respect to the T353 map and then calculate the corresponding 
HROs and ^ as a function of logio(A^H/cm“^). The intrinsic sta¬ 
tistical properties of the two maps are unchanged because the 
two maps are unchanged, only shifted. 

If the trend in ^ as a function of logio(A^H/cm“^) were aris¬ 
ing randomly, these trends would be unchanged even for signif¬ 
icant shifts. Instead the results of this experiment, illustrated in 
Fig. C.4 for the Taurus region, show that the trends tend to dis¬ 
appear with increasing values of the size of the shift. For shifts 
of about 1 °, the correlation between the magnetic field and the 
matter is still present, as expected from the results presented in 
Fig. 10, but for larger shifts the correlation is lost and the trend 
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Fig. C.3. Like Fig. C. 1 but for values of ^ (green) obtained using 
1000 random realizations of Q and U maps corresponding to a 
power spectrum P{k) oc k^^, with = -1.5, -2.5, and -3.5. 


does not survive. Over the many MCs studied, the nature of the 
trends at large shifts appears to be random. 

C.3. Projection effects 

We evaluate the statistical significance of the relative orientation 
between and VT 353 by considering the distribution of relative 
orientations between two vectors in 3D space compared to the 
distribution of relative orientation between their projections in 
2D. 

The projection of a vector v onto a plane normal to the unit 
vector h can be written as 

u = V - (v h)h. (C.l) 

Consider two unit vectors in 3D, vi and V 2 , separated by an 
angle a, such that 

cos a = Vi • V 2 . (C.2) 

The angle p between the projections of these two vectors onto a 
plane normal to w, Mi, and M 2 , can be written as 


COSyS = 


Ml • M2 

|Mi||M2| 

(vi •V2-(vi •n)(v 2 -n)) 


( Vi • Vi - (Vi • n)2 (V2 • V2 - (V2 • HY 


(C.3) 



Fig. C.4. Histogram shape parameter, as a function of 
logio(A^H/cm“^) calculated from the T353 map at 10 ' resolution 
in the Taurus region and the Stokes Q and U maps shifted in 
Galactic longitude and latitude by the values indicated, for A/ 
and A/?, respectively. 


Given a particular distribution of angles between the vectors Vi 
and V 2 , this expression, which is solved numerically, is useful for 
evaluating the resulting distribution of angles between the pro¬ 
jected vectors mi and M 2 . Without any loss of generality, we can 
assume that Vi is oriented along the axis of a spherical coordinate 
system, such that vi = k, and that V 2 is oriented at polar angle 6 
and azimuth 0 , such that V 2 = cos 0 sin 6i-\- sin 0 sin 6j-\- cos 6 k. 
Then, we can simulate a distribution of cos a by simulating a dis¬ 
tribution of cos 6 and generating 0 with a uniform distribution. 

We thus generate a set of vectors V 2 that follow a particu¬ 
lar distribution of cos a. We choose three examples: a uniform 
distribution of relative orientation between V 2 and Vi, V 2 vectors 
that are mostly parallel to Vi, and V 2 vectors that are mostly per¬ 
pendicular to Vi. The last two are Gaussian distributions centred 
at cos Of = 0 (for mostly parallel) or cos a = +1 (for mostly 
perpendicular, given the periodicity of a) with a dispersion 
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These distributions are shown in the top panel of Fig. C.5 for 
three values of 



jS [deg] 


Fig. C.5. Normalized distributions of relative orientations of ran¬ 
dom vectors in 3D (top) and of their projections in 2D (bottom). 
The curves correspond to mostly perpendicular vectors (red), 
uniform distribution of relative orientations (black), and mostly 
parallel vectors (green). The solid, dashed, and dotted lines cor¬ 
respond to dispersions of the relative orientation angle = 18°, 
36°, and 72°, respectively. 

Using Eq. (C.3) we calculate the distribution of the angles 
between the projected vectors and the results are shown in the 
bottom panel of Fig. C.5. Projection generally tends to make 
vector pairs look more parallel in 2D than they are in 3D, regard¬ 
less of the orientation in 3D, in agreement with the results of the 
Gaussian models presented in Planck Collaboration Int. XXXII 
(2014). Nevertheless, the distribution of relative orientations in 
2D is qualitatively similar (though not identical) to the distribu¬ 
tion in 3D. In particular, the signature of perpendicular orienta¬ 
tion is not erased. We conclude that the observation in 2D of 
close to perpendicular to the column density structures is a direct 
indicator of the perpendicular configuration of B with respect to 
the density structures in 3D. In contrast, the parallel alignment 
of B^ with the structure, or no signs of preferential orientation, 
while suggestive of parallel orientation in 3D, does not unam¬ 
biguously rule out the presence of some close-to-perpendicular 
orientations within the distribution in 3D. 


Appendix D: Alternative estimates of magnetic field 
strength 

Given the historic importance of the DCF method (Davis 1951; 
Chandrasekhar & Fermi 1953) and the related DCF-rSF method 


(Hildebrand et al. 2009), we have used them to estimate the 
magnetic field strength in each region. We also evaluated the 
mass-to-flux ratio, which is important in investigating the sta¬ 
bility against gravitational collapse. We now provide a critical 
assessment of the applicability of these results. 

D.1. Davis-Chandrasekhar-Fermi method 

The DCF method estimates the strength of B^ in a region using 
the dispersion of the polarization angles Assuming that the 
magnetic field is frozen into the gas and that the dispersion of the 
B^ orientation angles (or equivalently ^^) is due to transverse 
incompressible Alfven waves, then 

V4^^, (D.l) 

where is the dispersion of the radial velocity of the gas 
(Appendix D.3) and p the gas mass density. 

We calculate directly from Stokes Q and U using 

= ^<(A^)2) (D.2) 

and 

A-A= Iarctan(e(t/>-(e>t/, Q{Q) + {U)U) , (D.3) 

where (...) denotes an average over the selected pixels in each 
map (Planck Collaboration Int. XIX 2014). 

D.2. Davis-Chandrasekhar-Fermi plus structure function 
method 

As described by Hildebrand et al. (2009), the DCF-hSF method 
characterizes the magnetic field dispersion about local struc¬ 
tured fields by considering the difference in angle, = 

i//(x) - i//(x + i), between pairs of B^ vectors separated by dis¬ 
placements i in the plane of the sky. Assuming that the angle 
differences are statistically isotropic (i.e., they depend only on 
^ = \i\ and not on the orientation of i), they can binned by 
distance, £. From the N(^) pairs of B^ vectors for that bin, the 
square of the second-order structure function is 

/ 1 \ 

Sl({) = ) > (D-4) 

as introduced by Kobulnicky et al. (1994) and Falceta- 
Gon 9 alves et al. (2008). In terms of the Stokes parameters, each 
term in the sum can be written 

^ arctan {QtU^ - Q^Ui, QiQ^ + UiU^) , (D.5) 

where the subscripts v and i represent the central and displaced 
positions, respectively. 

Hildebrand et al. (2009) assume that B{x) is composed 
of a large-scale structured field, ^o(-^), and a random compo¬ 
nent, B^{x), which are statistically independent of each other. 
Assuming that ^o(-^) is a smoothly varying quantity, its contri¬ 
bution to S\{i) should increase in proportion to for distances 
that are much smaller than the scales at which Bq itself fluc¬ 
tuates. Additionally, assuming that turbulence occurs on scales 
that are small enough to completely decorrelate B^ for the range 
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Table D.l. Magnetic properties of the selected regions.^ 


Region 

U'vii 

[kms"^] 

[deg] 

[deg] 

b 

[rad] 

5DCF 

[bO] 

g DCF+SF 

(v«G] 

iDCF 

^obs 

1 DCF+SF 
^obs 

Taurus . 

1.2+0.5 

43+0.1 

23+0.05 

0.39+0.01 

13+5 

32+13 

0.4+0.4 

0 .2+0.1 

Ophiuchus . 

0.9+0.4 

29+0.3 

20+0.04 

0.36+0.01 

13+6 

25+11 

0.4+0.4 

0 .2+0.2 

Lupus . 

1.5+0.6 

46+0.7 

30+0.06 

0.52+0.01 

14+5 

29+11 

0.3+0.2 

0 .2+0.1 

Chamaeleon-Musca . 

1.0+0.4 

36+0.3 

23+0.05 

0.40+0.01 

12+5 

27+11 

0.4+0.3 

0 .2+0.2 

Corona Australis (CrA) 

0 .6+0.2 

59+0.1 

30+0.07 

0.52+0.01 

5+2 

12+ 5 

0.9+0.9 

0.3+0.3 

AquilaRift. 

1.9+0.6 

43+0.5 

23+0.09 

0.40+0.01 

20+6 

50+15 

0.3+0.2 

O.l+O.l 

Perseus. 

1.5+0.6 

38+0.3 

29+0.11 

0.50+0.01 

17+7 

30+11 

0.3+0.3 

0 .2+0.2 

IC5146 . 

1.7+0.6 

69+0.1 

49+0.11 

0.85+0.01 

11+4 

18+ 6 

0.5+0.3 

0.3+0.2 

Cepheus . 

1 .6+0.6 

43+0.2 

20+0.04 

0.35+0.01 

16+6 

47+18 

0.3+0.1 

O.l+O.O 

Orion. 

1.7+0.6 

36+0.1 

26+0.06 

0.45+0.01 

20+7 

38+14 

0.3+0.3 

0 .2+0.2 


^ Tabulated values are: the velocity dispersion, the dispersion of the polarization orientation angle, gf, the turbulent contribution to the 
angular dispersion, b, calculated from a linear fit to Eq (D.6); the magnetic field strengths, and calculated with the DCF method, 

of Eq. (D.l), and the DCF+SF method, of Eq. (D.IO), with b values obtained from a fit to Eq. (D.6) in the range 50'< £ < 200'; and their 
corresponding observed mass-to-flux ratios, and The reported uncertainties are from the appropriate propagation of errors. 


of scales probed by the displacements £, Hildebrand et al. (2009) 
derived the approximation 

slit) = b^ + , (D.6) 


where the two terms on the right-hand side give the contribu¬ 
tions from the random and large-scale magnetic fields, respec¬ 
tively. These assumptions are not necessarily valid for the range 
of scales in the Planck data. 

Hildebrand et al. (2009) used Eq. (D.6) to estimate b^ as the 
intercept and then related b to the ratio of the random to the 
large-scale magnetic field strength, both projected onto the plane 
of the sky, through 





b 

V2 - ’ 


(D.7) 


/ 9 \l/2 

where stands for the root mean square (rms) variations 

about the large-scale magnetic field, 5o,±- The same assumptions 
that result in Eq. (D.l) - i.e., considering only incompressible 
and isotropic turbulence, magnetic fields frozen into the gas, and 
dispersion of the orientation originating in transverse incom¬ 
pressible Alfven waves - lead to 





(D.8) 


where 

Va,x = (D.9) 

is the speed of the transverse incompressible Alfven waves. 
Combining Eqs. (D.7)-(D.9) results in 


dDCF+SF _ dDCF+SF 
- ^o,± 


^jAnp 


O-v,, V 2-^2 

b 


(D.IO) 


D.3. Calculation 

The estimates of velocity dispersion (Tv|| and mass density p are 
the same in the two methods (DCE and DCE-rSE). 

We obtain from the most complete CO emission-line sur¬ 
vey of the Milky Way currently available, that of Dame et al. 


(2001). This data set consists of 488 000 spectra that beam- 
sample (1/8°) a set of MCs and the Galactic plane over a strip 
from 4° to 10° wide in latitude. We find that material in the MCs 
has centroid velocities in the range -10 < V||/(kms“^) < 10 
in the first moment map. Eurthermore, for the calculations be¬ 
low, we need to select pixels with sufficient S/N in the polar¬ 
ization observations and so use the second criterion described 
in Appendix A.2. This combination of cuts selects the areas il¬ 
lustrated in Eig. D.l. Over these areas we calculate the average 
of the velocity dispersion from the second moment map and use 
this value as . This is tabulated along with other properties in 
Table D.l. 

The mass density, p = yu^mp, is the product of the proton 
mass, mp, the mean molecular weight per hydrogen molecule, 
p = 2.8, and the mean number density, n. We require a value of n, 
which for the discussions here we approximate to be 100 cm“^. 
This is a typical value for MCs (Draine 2011) and in rough 
agreement with the column densities and cloud sizes presented 
in Table 1 . In practice, n varies from one cloud to the next, and its 
estimation involves the assumption of a particular cloud geom¬ 
etry, resulting in additional uncertainties that are not considered 
in this study. 

Eor the DCE method, the dispersion of the orientation angle 
corresponds to the centred second moment of the if/ distribu¬ 
tion evaluated in pixels within the selected region. Eollowing the 
analysis presented in Planck Collaboration Int. XIX (2014) for 
S{£), the variance on can be expressed as 



where if/i and cr^|^. are the orientation angle and uncertainty for 
each polarization pseudo-vector, N is the number of if/ measure¬ 
ments, gi = ij/i - (i//), and {i//) is the average orientation angle in 
each region. Like other quadratic functions, cr'^^ is biased posi¬ 
tively when noise is present, leading to an overestimation of this 
quantity. However, given that we limit our analysis to polariza¬ 
tion measurements with high S/N and that the uncertainties in 
(Tv|| are considerably larger, the bias correction does not have a 
significant effect on our estimates. We apply Eq. (D.l) to esti¬ 
mate and propagate the errors to obtain the values listed in 
Table D.l. 
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[km/s] ^ [km/s] ^ [km/s] 



180 176 172 168 164 0.98 -2.51 -6.00 -9.49 -12.98 346 343 340 337 334 

L [deg] I [deg] L [deg] 



308.5 304.3 300.0 295.7 291.5 

Z [deg] 


6.70 3.35 0.00 -3.35 -6.70 

I [deg] 



33 


30 27 24 

Z [deg] 


21 163.9 161.4 159.0 156.6 154.1 

Z [deg] 



[km/s] Q [km/s] 



118.5 114.3 110.0 105.7 101.5 220.6 216.3 212.0 207.7 203.4 

Z [deg] Z [deg] 


Fig.D.l. Line-of-sight centroid velocity vy maps inferred from CO emission (Dame et al. 2001). Areas shown have sufficient S/N 
in the polarization observations (second polarization criterion in Appendix A.2) for the DCF and DCF+SF analyses and lie in the 
range -10 < V||/(kms“^) < 10. 
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I [orcmin] 



0 50 100 150 200 

I [orcmin] 


Fig.D.2. Structure function S 2 {^) calculated from the Q and U maps of the selected regions following Eq. (D.4). The black lines 
indicate the fits given by Eq. (D.6). The vertical dashed line marks the common 10' resolution of the data used in the analysis; the 
data are correlated for low values of i causing the drop in 5 2 (^). 


Eor the DCE+SE method, to estimate the parameter b, we 
first calculate for each of the regions. The resolution of 
the data used is 10'. We evaluate Sl(£) in steps of 3.'44 for lags 
0' < ^ < 34.'4 and in steps of 34.'4 for lags 40' < ^ < 200'. Eor 
each considered lag is evaluated pixel by pixel, consid¬ 

ering all the pixels located in an annulus with radius £. The term 
Ai//xj is only considered in the calculation of S if there are at 
least three pixels in the annulus. As anticipated, the model from 
Eq. (D.6) does not agree with the data on all sampled scales, so 
the range of scales considered is limited to £ above the resolu¬ 


tion of the data and between 50' and 200', where the behaviour 
of is approximately linear in We then make a linear 
fit of Sl(^) using as the independent variable and estimate 
the value of b and its uncertainty, cr^. The results are plotted in 
Eig. D.2. 

Using the calculated b values,^ we apply Eq. (D.IO) to es¬ 
timate and the results are listed in Table D.l. As ex- 


^ For Eq. (D.IO) b is required to be in radians, and so those are also 
listed in Table D.l. 
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pected, the DCF+SF method produces estimates of that are 
about twice as large as those obtained with the DCF method. 
However, it is worth noting that because of the discrepancy be¬ 
tween the model and the observations, the value of b depends on 
the selected ^-range and so propagates into a different value of 

^ DCF+SF 


DA. Mass-to-flux ratios 

The critical value for the mass that can be supported against 
gravity by a magnetic flux O can be estimated to first or¬ 
der for a uniform disk from (M/0)crit = l/(27rG^/^) (Nakano 
& Nakamura 1978). The precise value of the right-hand side 
changes for different geometries (e.g., Spitzer 1968; McKee 
et al. 1993). Stability can be assessed using 


(M/O) 

(M/0)erit 


= 7.6 X 10“2' 




(D.12) 


where and are the H 2 column density and magnetic 
field strength along a magnetic flux tube (Crutcher et al. 2004). A 
cloud is supercritical and prone to collapse under its own gravity, 
when T > 1; otherwise, when T < 1, the cloud is sub-critical, 
magnetically supported against gravitational collapse. 

What is observable is Tobs, in which is replaced by 

We evaluate /fobs by combining the value of (A^H 2 ) com¬ 
puted from the integrated CO line emission and the conversion 
factor Xco = (1.8 + 0.3) x 10^^ cm“^ K“^ km“^ s (Dame et al. 
2001) and estimated with the DCF and DCF-i-SF methods. 
The Xco factor may show cloud-to-cloud or regional variations 
(Draine 2011), but we consider that these are not significant in 
comparison to the uncertainties involved in the estimation of Bj^. 
The calculated values of Tobs are listed in Table D.l. They are 
consistent with being less than unity. 

From this we might obtain A by judicious deprojection, since 


A = (N^^/Nh,) X dobs = /dp dobs. (D.13) 


Here, B^ is always less than B^^\ pushing /dp below unity. The 
situation for column density depends on the geometry of the 
structure relative to the magnetic field. For a structure with the 
magnetic field along the short axis, < Ah 2 , again lowering 

/dp- 

Statistically, the mean mass-to-flux ratio can be related to the 
observed value by assuming a particular geometry of the cloud, 
an ellipsoid with equatorial radius a and centre-to-pole distance 
c, and a magnetic field oriented along the polar axis of the ellip¬ 
soid. For an oblate spheroid, flattened perpendicular to the ori¬ 
entation of the magnetic field, /dp = 1/3, yielding 

- McosB 1 

(D.14) 

Jo O/ sm/3 3 


where p is the inclination angle with respect to the line of sight, 
and the sin p dependence in the flux comes from B = B^ sin p, 
(Crutcher et al. 2004). For a sphere there is no cosyS dependence 
and /dp = 1/2. For a prolate spheroid elongated along the orien¬ 
tation of the field, the mass is multiplied by sinyS instead of cos p, 
resulting in /dp = 3/4. Investigating which geometry, if any, is 
most relevant to the actual MCs and the magnetized structures 
detected within them is obviously important before firm conclu¬ 
sions can be drawn regarding gravitational instability. 


D.5. Discussion 

Our estimates of the magnetic field strengths in the MCs 
analysed and the mass-to-flux ratios, presented in Table D.l, 
stem from the classic calculation presented by Chandrasekhar 
& Fermi (1953) and an updated interpretation presented by 
Hildebrand et al. (2009). Both of these methods assume very 
specific conditions, so (as we discuss below) this limits the con¬ 
clusions that can be drawn from these estimates of the magnetic 
field strength. The deduced mass-to-flux ratios suggest that the 
clouds are potentially magnetically sub-critical, but we again 
need to be cautious about drawing conclusions from this appli¬ 
cation of the DCF and DCF-rSF analyses alone. 

In the case of the DCF method, the values of B^^^ are ob¬ 
tained by assuming that the structure of the magnetic field is 
the product only of incompressible Alfven waves, where the dis¬ 
placements are perpendicular to the direction of propagation. 
This is not the case for turbulence in MCs where the random 
component of the magnetic field can have any orientation. The 
dispersion measured about mean fields, assumed to be straight 
in DCF, may be much larger than should be attributed to MHD 
waves or turbulence, leading to an overestimation of and to 
low values of Bj^. Moreover true interstellar turbulence in MCs 
involves not only incompressible Alfven waves, but also com¬ 
pressible magneto-sonic waves, which do not satisfy Eq. D.l. 
Furthermore, depending on the scales examined, the magnetic 
field may have structures due to effects such as differential rota¬ 
tion, gravitational collapse, or expanding Hii regions. 

In the case of the DCF-i-SF method, the values of ^^cf+sf 
are obtained by assuming a very specific model of the magnetic 
field, which is in principle just a first-order approximation. First, 
this model assumes that the effect of the large-scale structured 
magnetic field, Bq, is to cause the square of the second-order 
structure function, 5 2 (^), to increase as This corresponds to a 
very specific correlation function for Bq. Second, this model as¬ 
sumes that the dispersion of the random component of the field. 
By, is scale-independent, which is not realistic for the range of 
scales probed by Planck (Elmegreen & Scalo 2004; Hennebelle 
& Falgarone 2012). 

In addition, the magnetic field orientation deduced from the 
polarization angle in a particular direction is not generally that of 
the field at a single point along the line of sight. The observed po¬ 
larization is the average of various field pseudo-vectors weighted 
by local dust emission along the line of sight. The net effect of 
the integration of multiple uncorrelated components along the 
line of sight is an observed dispersion of the polarization angle 
that is smaller than the true 3D dispersion of the magnetic field 
orientation, thus leading to an overestimation of B^. Myers & 
Goodman (1991) presented an analysis of this effect in terms of 
the number of correlation lengths of the magnetic field along the 
line of sight through a cloud, which they calculated empirically. 

Houde et al. (2009) presented an extension of the DCF-rSF 
method that includes the effect of signal integration along the 
line of sight and across the area subtended by the telescope 
beam. The extended method, also implemented in Houde et al. 
(2011, 2013), is based on the identification of the magnetized 
turbulence correlation length (^) by means of the structure func¬ 
tion of the polarization angles. In the case of the Planck 353 GHz 
observations, the angular resolution is not sufficient to identify 6 
and the corrections would have to rely on rough estimates of this 
value and the depth of integration (A). Following equation 29 in 
Houde et al. (2009), rough estimates 6 ^ 0.2 pc and A 10 pc 
result in a few correlation lengths across the beam, correspond¬ 
ing to correction factors around 0.4. Coincidentally, such cor- 
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rection factors lead to values of ^ dcf+sf those of 

in Table D.l. But note that this correction relies on specific as¬ 
sumptions on the nature of the turbulence correlation function, it 
does not circumvent the necessity to observe the dust polarized 
emission with higher angular resolution to fully characterize the 
magnetic field. 

MHD simulations provide a potentially useful guide to what 
modifications might be introduced into the DCF formula to al¬ 
low for inhomogeneity and line-of-sight averaging. Using syn¬ 
thetic observations of MHD simulations, Ostriker et al. (2001) 
showed that correcting Eq. (D.l) by a factor C ~ 0.5 provides a 
good approximation to the actual magnetic field strength in cases 
where < 25°. However, the effect of nonlinear amplitudes 
is uncertain (Zweibel 1996), and the method fails for values of 

> 25°, which is the case for all of the regions in this study 
(Table D.l). Falceta-Gon 9 alves et al. (2008) propose a method 
that is potentially valid for any value of 9 *^, based on a fit to the 
^DCF yaiues obtained from maps at different resolutions, again 
concluding that the field should be lower than estimated with 
Eq. (D.l). However, this approach was not tested in MHD simu¬ 
lations that include gravity, which is the critical process that we 
aim to evaluate by using the DCF method. 

Another strong assumption is that the behaviour of the veloc¬ 
ity and the magnetic field are represented well by the observed 
quantities and 9 *^ for a particular set of scales, which might 
not necessarily be the case. Even if the power spectra of v and 
B are comparable in 3D, the integration along the line of sight 
is different for the two quantities. The dispersion is based on 
the emission-line profile V|| tracing a gas species, and while the 
emission is directly integrated along the line of sight, the line 
profile is possibly affected by radiative transfer and excitation 
effects. The tracer of the magnetic field is the optically-thin po¬ 
larized submillimetre emission of dust, and both the polarization 
and are projected and integrated along the line of sight as 
pseudo-vectors. 

Finally, in the estimates of and a common 

mean density p has been adopted, while in practice p is different 
from cloud to cloud. Direct estimation of the values of p from 
the measured column densities (Ah) and (AH 2 ) relies heavily on 
the geometrical modelling of the cloud and the filling factors of 
each species, introducing uncertainties that affect the calculated 
values of B^. 

Given the line-of-sight integration and the fact that > 25° 
uniformly, the calculated values could be considered as upper 
limits to the actual B^. However, other shortcomings, such as 
the assumptions about the correlation structure and the uncer¬ 
tainties in the determination of the density, do not necessarily 
bias the estimate of the magnetic field strength towards high val¬ 
ues. In conclusion, the values presented in Table D.l should be 
viewed only as a reference and only applied with caution, given 
the many assumptions in both methods at the scales considered. 
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